---
title: "Mosques and Religious Attitudes"
author: |
  | Nathanael Gratias Sumaktoyo
  | National University of Singapore
  | nathanael.sumaktoyo@nus.edu.sg
output:
  html_document:
    toc: yes
    toc_float: yes
    toc_collapsed: yes
    toc_depth: 3
    number_sections: yes
    theme: lumen
    df_print: paged
---

version: `r format(Sys.time(), '%d %B, %Y -- %H:%M:%S')`

```{r, echo=FALSE, warning=FALSE, results='hide', include=FALSE}
rm(list=ls())
source("extras/packages.R")
source("extras/functions.R")
knitr::opts_chunk$set(rows.print=100)

```


# Figure-2 


```{r, echo=FALSE, warning=FALSE}

### mosque data
load("data/mosque_data.RData")

byyear = mosque %>%
  group_by(year) %>%
  dplyr::summarise(n = n ())

byyear = subset(byyear, year %in% 1945:2022)

sums = cumsum(byyear$n)
sums = data.frame(sums)
sums$year = byyear$year
sums = subset(sums, year %in% 1990:2022)
sums$sums = sums$sums / 1000

p = ggplot(sums, aes(x = year, y = sums)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  geom_vline(xintercept = 1998, lty = 2) +
  geom_label(aes(x = 1998, y = 120), label = "Collapse of\nSuharto's Regime") +
  geom_label(aes(x = 1998, y = 70), label = "Start of\nDemocratization") +
  geom_vline(xintercept = 2014, lty = 2) +
  geom_label(aes(x = 2014, y = 120), label = "SIMAS\nIntroduced") +
  scale_x_continuous("Year",
                     breaks = seq(1990,2022,2)) +
  scale_y_continuous("Total Number of Mosques\n(in Thousands)", 
                     breaks = seq(0, 280, 20),
                       limits = c(0, 280)) + 
  labs(title = "Number of Mosques Registered on SIMAS (1990 - 2022), in Thousands")

p

ggsave(p, filename = "Figure-2.pdf")
```


```{r, echo=FALSE, message=FALSE, warning=FALSE, results='hide'}

# https://easystats.github.io/parameters/articles/model_parameters_robust.html#cluster-robust-covariance-matrix-estimation-sandwich

# https://strengejacke.github.io/ggeffects/articles/practical_robustestimation.html#predictions-with-cluster-robust-standard-errors

rm(list = setdiff(ls(), lsf.str()))

source("extras/packages.R")
source("extras/functions.R")

### kecamatan data
load("data/dat_kec_long.RData")


dat = subset(dat, subset = (religion_07 == 1) & (religion_14 == 1) & 
                            (census10_prop_muslim >= .5) & (mig_status == "stayer"))

dat$year = car::recode(dat$year, "'07' = '2007'; '14' = '2014'")
dat$year = factor(dat$year, levels = c("2007", "2014"))

dat$shalatlab3pt = car::recode(dat$shalat3pt,
                               "1 = 'less than five';
                               2 = 'five';
                               3 = 'more than five'")
dat$shalatlab3pt = factor(dat$shalatlab3pt, 
                          levels = c("five", "less than five", "more than five"))




```

Kecamatans and respondents by treatment status:

```{r, echo=FALSE}
dat$treat = ifelse(dat$change_mo_07_13 > 0, "New Mosque(s)", "No New Mosque")
dat$treat = factor(dat$treat, levels = c("No New Mosque", "New Mosque(s)"))

table(dat$treat) / 2
prop.table(table(dat$treat)) 
```



# Figure-3

```{r, echo=FALSE, fig.width=9}

coefs = read.csv("./coefs_robust_placebo/coefs_main_models.csv", header=TRUE)

coefs$low99 = coefs$coef - qnorm(.995) * coefs$se
coefs$hi99 = coefs$coef + qnorm(.995) * coefs$se

coefs$low95 = coefs$coef - qnorm(.975) * coefs$se
coefs$hi95 = coefs$coef + qnorm(.975) * coefs$se

coefs$low90 = coefs$coef - qnorm(.95) * coefs$se
coefs$hi90 = coefs$coef + qnorm(.95) * coefs$se


dvs = c("Object: Live in Same Village (a)", 
        "Object: Live in Same Neighborhood (a)", 
        "Object: Live in Same House (a)",
        "Object: Worship House (b)", 
        "Object: Intermarriage (a)", 
        "Prefer Muslim Candidate (c)", 
        "Prefer Co-Ethnic Candidate (c)",
        "Trust Fellow Muslims (d)", 
        "Trust Co-Ethnics (a)", 
        "Helping Neighbors (a)")

coefs$category = c(rep("Outgroup\nRejection", 5),
                   rep("Political\nPreferences", 2),
                   rep("Ingroup\nAttitudes", 3))
  
coefs$category = factor(coefs$category,
                       levels = c("Outgroup\nRejection", "Political\nPreferences", 
                                  "Ingroup\nAttitudes"))



coefs$dv = dvs
coefs$dv = factor(coefs$dv,
                  levels = coefs$dv[length(dvs):1])
idx = which(coefs$p < (.05/nrow(coefs)))

p <- ggplot(coefs, 
       aes(x = dv, y = coef)) +
        geom_hline(yintercept = 0, 
                   colour = "red", lty = 2.5) +
        geom_point(aes(x = dv, 
                    y = coef), cex=2) + 
        geom_linerange(aes(x = dv, 
                     ymin = low95,
                     ymax = hi95),
                   linewidth = 1) +
        # geom_linerange(aes(x = dv, 
        #              ymin = low99,
        #              ymax = hi99),
        #            linewidth = .5) + 
        scale_x_discrete(name = "Dependent Variables") +
        scale_y_continuous(name = "ATT") + 
        labs(title = "Effects of New Mosques across Dependent Variables",
             subtitle = "(95% Confidence Interval)",
             caption = "a: 16,238 respondents in 1,080 districts     c: 16,144 respondents in 1,077 districts\nb: 16,237 respondents in 1,080 districts     d: 16,236 respondents in 1,080 districts\n") +
        coord_flip() +
        facet_grid(rows = vars(category),
             scales = "free", space = "free_y") +
        theme_bw() +
        theme(axis.text=element_text(size=11))

p

ggsave(filename = "Figure-3.pdf", p, width = 10, height = 8)

```

# Figure-4

```{r, echo=FALSE, warning=FALSE, fig.width=9}
spec0 = read.csv("./coefs_robust_placebo/coefs_main_models.csv", header=TRUE)
spec1 = read.csv("./coefs_robust_placebo/coefs_placebo_inspace.csv", header=TRUE)
spec2 = read.csv("./coefs_robust_placebo/coefs_placebo_intime.csv", header=TRUE)
spec3 = read.csv("./coefs_robust_placebo/coefs_placebo_movers.csv", header=TRUE)


df = rbind(spec0[,1:4], spec1, spec2, spec3)
  

dvs = c("Object: Live in Same Village", 
        "Object: Live in Same Neighborhood", 
        "Object: Live in Same House",
        "Object: Worship House", "Object: Intermarriage", 
        "Prefer Muslim Candidate", "Prefer Co-Ethnic Candidate",
        "Trust Fellow Muslims", "Trust Co-Ethnics", 
        "Helping Neighbors")
df$dv = dvs
df$dv = factor(df$dv, levels = dvs[length(dvs):1])


mods = c("Main Models",
         "In-Space Placebo",
         "In-Time Placebo",
         "Analysis of Movers")
ndv = length(unique(df$dv))
df$modellab = rep(mods, times = rep(ndv, length(mods)))
     
df$modellab = factor(df$modellab, levels = mods)

df$category = NA
df$category[df$dv %in% c("Object: Live in Same Village", 
                          "Object: Live in Same Neighborhood", 
                          "Object: Live in Same House",
                          "Object: Worship House", 
                         "Object: Intermarriage")] = "Outgroup\nRejection"

df$category[df$dv %in% c("Prefer Muslim Candidate", 
                         "Prefer Co-Ethnic Candidate")] = "Political\nPreferences"

df$category[df$dv %in% c("Trust Fellow Muslims", 
                         "Trust Co-Ethnics", 
                         "Helping Neighbors")] = "Ingroup\nAttitudes"


df$category = factor(df$category,
                       levels = c("Outgroup\nRejection", "Political\nPreferences", 
                                  "Ingroup\nAttitudes"))

df$low95 = df$coef - qnorm(.975)*df$se
df$hi95 = df$coef + qnorm(.975)*df$se

df$col = "darkgrey"
df$col[(df$coef < 0) & (df$p<=.05)] = "red"
df$col[(df$coef > 0) & (df$p<=.05)] = "blue"

df$alph = .6
df$alph[df$p<=.05] = 1



d = subset(df, modellab == "Main Models")
xmin = min(d$low95)
xmax = max(d$hi95)
p1 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "pt"))


d = subset(df, modellab == "In-Space Placebo")
xmin = min(d$low95)
xmax = max(d$hi95)
p2 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == "In-Time Placebo")
xmin = min(d$low95)
xmax = max(d$hi95)
p3 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == "Analysis of Movers")
xmin = min(d$low95)
xmax = max(d$hi95)
p4 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(plot.margin = margin(0, 0, 0, 2, "pt")) 


p <- p1 + p2 + p3 + p4 + plot_layout(nrow=1)

p
ggsave(p, filename = "Figure-4.pdf")


```



# Figure-5


```{r, echo=FALSE, warning=FALSE, fig.width=9}
spec0 = read.csv("./coefs_robust_placebo/coefs_main_models.csv", header=TRUE)
spec1 = read.csv("./coefs_robust_placebo/coefs_specification_1.csv", header=TRUE)
spec2 = read.csv("./coefs_robust_placebo/coefs_specification_2.csv", header=TRUE)
spec3 = read.csv("./coefs_robust_placebo/coefs_specification_3.csv", header=TRUE)
spec4 = read.csv("./coefs_robust_placebo/coefs_specification_4.csv", header=TRUE)
spec5 = read.csv("./coefs_robust_placebo/coefs_specification_5.csv", header=TRUE)
spec6 = read.csv("./coefs_robust_placebo/coefs_specification_6.csv", header=TRUE)


df = rbind(spec0[,1:4], spec1, spec2, spec3, spec4, spec5, spec6)
  

dvs = c("Object: Live in Same Village", 
        "Object: Live in Same Neighborhood", 
        "Object: Live in Same House",
        "Object: Worship House", 
        "Object: Intermarriage", 
        "Prefer Muslim Candidate", 
        "Prefer Co-Ethnic Candidate",
        "Trust Fellow Muslims", 
        "Trust Co-Ethnics", 
        "Helping Neighbors")
df$dv = dvs
df$dv = factor(df$dv, levels = dvs[length(dvs):1])


mods = c("Binary DV\nNo Covariates\n(Main Models)",
         "Binary DV\nwith Covariates",
         "Ordinal DV\nNo Covariates",
         "Ordinal DV\nwith Covariates",
         "Main Models,\nDistricts (n > 10)",
         "Main Models,\nZero or One\nNew Mosque",
         "Main Models\nwith Matching")
ndv = length(unique(df$dv))
df$modellab = rep(mods, times = rep(ndv, length(mods)))
     
df$modellab = factor(df$modellab, levels = mods)

df$category = NA
df$category[df$dv %in% c("Object: Live in Same Village", 
                          "Object: Live in Same Neighborhood", 
                          "Object: Live in Same House",
                          "Object: Worship House", 
                         "Object: Intermarriage")] = "Outgroup\nRejection"

df$category[df$dv %in% c("Prefer Muslim Candidate", 
                         "Prefer Co-Ethnic Candidate")] = 
  "Political\nPreferences"

df$category[df$dv %in% c("Trust Fellow Muslims", 
                         "Trust Co-Ethnics", 
                         "Helping Neighbors")] = "Ingroup\nAttitudes"

df$category = factor(df$category,
                       levels = c("Outgroup\nRejection",
                                  "Political\nPreferences", 
                                  "Ingroup\nAttitudes"))

df$low95 = df$coef - qnorm(.975)*df$se
df$hi95 = df$coef + qnorm(.975)*df$se

df$col = "darkgrey"
df$col[(df$coef < 0) & (df$p<=.05)] = "red"
df$col[(df$coef > 0) & (df$p<=.05)] = "blue"

df$alph = .6
df$alph[df$p<=.05] = 1



d = subset(df, modellab == mods[1])
xmin = min(d$low95)
xmax = max(d$hi95)
p1 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "pt"))


d = subset(df, modellab == mods[2])
xmin = min(d$low95)
xmax = max(d$hi95)
p2 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[3])
xmin = min(d$low95)
xmax = max(d$hi95)
p3 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[4])
xmin = min(d$low95)
xmax = max(d$hi95)
p4 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[5])
xmin = min(d$low95)
xmax = max(d$hi95)
p5 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 

d = subset(df, modellab == mods[6])
xmin = min(d$low95)
xmax = max(d$hi95)
p6 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[7])
xmin = min(d$low95)
xmax = max(d$hi95)
p7 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(plot.margin = margin(0, 0, 0, 2, "pt")) 


p <- p1 + p2 + p3 + p4 + p5 + p6 + p7 + plot_layout(nrow=1)

p
ggsave(p, filename = "Figure-5.pdf", width = 10)

```




# Figure-6

```{r, echo=FALSE}

# write.csv(coefs, file = "./coefs_robust_placebo/coefs_religiosity.csv", row.names = FALSE)
coefs = read.csv("./coefs_robust_placebo/coefs_religiosity.csv", header=TRUE)

coefs$low99 = coefs$coef - qnorm(.995) * coefs$se
coefs$hi99 = coefs$coef + qnorm(.995) * coefs$se

coefs$low95 = coefs$coef - qnorm(.975) * coefs$se
coefs$hi95 = coefs$coef + qnorm(.975) * coefs$se

coefs$low90 = coefs$coef - qnorm(.95) * coefs$se
coefs$hi90 = coefs$coef + qnorm(.95) * coefs$se


coefs$dv = factor(coefs$dv,
                  c("Shalat Five Times or More", 
                    "Self-Rated Religiosity")
                  )


p1 <- ggplot(coefs, 
       aes(x = dv, y = coef)) +
        geom_hline(yintercept = 0, 
                   colour = "red", lty = 2) +
        geom_point(aes(x = dv, 
                    y = coef), cex = 2.5) + 
        geom_linerange(aes(x = dv, 
                     ymin = low95,
                     ymax = hi95),
                   linewidth = 1) +
        # geom_linerange(aes(x = dv, 
        #              ymin = low99,
        #              ymax = hi99),
        #            linewidth = .5) + 
        scale_x_discrete(name = "Dependent Variables") +
        scale_y_continuous(name = "ATT") + 
        labs(title = "Effects of New Mosques across Religiosity Measures",
             subtitle = "(95% Confidence Interval)") +
        coord_flip() +
        theme_bw()

p1
```


```{r, echo=FALSE}

# write.csv(coefs, file = "./coefs_robust_placebo/coefs_mushallas.csv", row.names = FALSE)
coefs = read.csv("./coefs_robust_placebo/coefs_mushallas.csv", header=TRUE)

coefs$low99 = coefs$coef - qnorm(.995) * coefs$se
coefs$hi99 = coefs$coef + qnorm(.995) * coefs$se

coefs$low95 = coefs$coef - qnorm(.975) * coefs$se
coefs$hi95 = coefs$coef + qnorm(.975) * coefs$se

coefs$low90 = coefs$coef - qnorm(.95) * coefs$se
coefs$hi90 = coefs$coef + qnorm(.95) * coefs$se

coefs$dv = c("Object: Live in Same Village", 
        "Object: Live in Same Neighborhood", 
        "Object: Live in Same House",
        "Object: Worship House", "Object: Intermarriage", 
        "Prefer Muslim Candidate", "Prefer Co-Ethnic Candidate",
        "Trust Fellow Muslims", "Trust Co-Ethnics", 
        "Helping Neighbors")

coefs$dv = factor(coefs$dv,
                  c("Helping Neighbors", 
                    "Trust Co-Ethnics", 
                    "Trust Fellow Muslims", 
                    "Prefer Co-Ethnic Candidate",
                    "Prefer Muslim Candidate", 
                    "Object: Intermarriage",
                    "Object: Worship House",
                    "Object: Live in Same House",
                    "Object: Live in Same Neighborhood", 
                    "Object: Live in Same Village"
                    )
                  )

coefs$category = c(rep("Outgroup\nRejection", 5),
                   rep("Political\nPreferences", 2),
                   rep("Ingroup\nAttitudes", 3))
  
coefs$category = factor(coefs$category,
                       levels = c("Outgroup\nRejection", "Political\nPreferences", 
                                  "Ingroup\nAttitudes"))


p2 <- ggplot(coefs, 
       aes(x = dv, y = coef)) +
        geom_hline(yintercept = 0, 
                   colour = "red", lty = 2) +
        geom_point(aes(x = dv, 
                    y = coef), cex = 2.5) + 
        geom_linerange(aes(x = dv, 
                     ymin = low95,
                     ymax = hi95),
                   linewidth = 1) +
        # geom_linerange(aes(x = dv, 
        #              ymin = low99,
        #              ymax = hi99),
        #            linewidth = .5) + 
        scale_x_discrete(name = "Dependent Variables") +
        scale_y_continuous(name = "ATT") + 
        labs(title = "Effects of New Musallas across Dependent Variables",
             subtitle = "(95% Confidence Interval)") +
        coord_flip() +
        facet_grid(rows = vars(category),
             scales = "free", space = "free_y") +
        theme_bw()

p2

ggsave(p1, filename = "Figure-6a.pdf")
ggsave(p2, filename = "Figure-6b.pdf")

```


# Figure-7

```{r, echo=FALSE, warning=FALSE, fig.width=9}
spec1 = read.csv("./coefs_int_ind_cov/friday.csv", header=TRUE)
spec2 = read.csv("./coefs_int_ind_cov/numuha.csv", header=TRUE)
spec3 = read.csv("./coefs_int_kec_cov/propmuslim.csv", header=TRUE)
spec4 = read.csv("./coefs_int_kec_cov/propwakaf.csv", header=TRUE)
spec5 = read.csv("./coefs_int_kec_cov/nmo07.csv", header=TRUE)
spec6 = read.csv("./coefs_int_kec_cov/numuha.csv", header=TRUE)

df = rbind(spec1, spec2, spec3, spec4, spec5, spec6)


mods = c("Friday Interview\n(Individual)",
         "Muslim Organization\nAffiliation\n(Individual)",
         "Prop. Muslim\n(District)",
         "Prop. Waqf\nMosques\n(District)",
         "Log(N Mosques\nin 2007)\n(District)",
         "Proportion of\nMuslim Organization\nAffiliation\n(District)")


dvs = c("Object: Live in Same Village", 
        "Object: Live in Same Neighborhood", 
        "Object: Live in Same House",
        "Object: Worship House", "Object: Intermarriage", 
        "Prefer Muslim Candidate", "Prefer Co-Ethnic Candidate",
        "Trust Fellow Muslims", "Trust Co-Ethnics", 
        "Helping Neighbors")
df$dv = dvs
df$dv = factor(df$dv, levels = dvs[length(dvs):1])


ndv = length(unique(df$dv))
df$modellab = rep(mods, times = rep(ndv, length(mods)))
     
df$modellab = factor(df$modellab, levels = mods)

df$category = NA
df$category[df$dv %in% c("Object: Live in Same Village", 
                          "Object: Live in Same Neighborhood", 
                          "Object: Live in Same House",
                          "Object: Worship House", 
                         "Object: Intermarriage")] = "Outgroup\nRejection"

df$category[df$dv %in% c("Prefer Muslim Candidate", 
                         "Prefer Co-Ethnic Candidate")] = 
  "Political\nPreferences"

df$category[df$dv %in% c("Trust Fellow Muslims", 
                         "Trust Co-Ethnics", 
                         "Helping Neighbors")] = "Ingroup\nAttitudes"


df$category = factor(df$category,
                       levels = c("Outgroup\nRejection",
                                  "Political\nPreferences", 
                                  "Ingroup\nAttitudes"))


df$low95 = df$coef - qnorm(.975)*df$se
df$hi95 = df$coef + qnorm(.975)*df$se

df$col = "darkgrey"
df$col[(df$coef < 0) & (df$p<=.05)] = "red"
df$col[(df$coef > 0) & (df$p<=.05)] = "blue"

df$alph = .6
df$alph[df$p<=.05] = 1



d = subset(df, modellab == mods[1])
xmin = min(d$low95)
xmax = max(d$hi95)
p1 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "pt"))


d = subset(df, modellab == mods[2])
xmin = min(d$low95)
xmax = max(d$hi95)
p2 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[3])
xmin = min(d$low95)
xmax = max(d$hi95)
p3 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[4])
xmin = min(d$low95)
xmax = max(d$hi95)
p4 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == mods[5])
xmin = min(d$low95)
xmax = max(d$hi95)
p5 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 



d = subset(df, modellab == mods[6])
xmin = min(d$low95)
xmax = max(d$hi95)
p6 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(plot.margin = margin(0, 0, 0, 2, "pt")) 

p <- p1 + p2 + p3 + p4 + p5 + p6 + plot_layout(nrow=1)

p
ggsave(p, filename = "Figure-7.pdf", width = 10)

```



# Figure-8

```{r, echo=FALSE, warning=FALSE, fig.width=9}
spec1 = read.csv("./coefs_int_ind_cov/edu.csv", header=TRUE)
spec2 = read.csv("./coefs_int_ind_cov/cognitive.csv", header=TRUE)
spec3 = read.csv("./coefs_int_ind_cov/infosources.csv", header=TRUE)
spec4 = read.csv("./coefs_int_kec_cov/medianedu.csv", header=TRUE)


df = rbind(spec1, spec2, spec3, spec4)
  

dvs = c("Object: Live in Same Village", 
        "Object: Live in Same Neighborhood", 
        "Object: Live in Same House",
        "Object: Worship House", "Object: Intermarriage", 
        "Prefer Muslim Candidate", "Prefer Co-Ethnic Candidate",
        "Trust Fellow Muslims", "Trust Co-Ethnics", 
        "Helping Neighbors")
df$dv = dvs
df$dv = factor(df$dv, levels = dvs[length(dvs):1])


mods = c("Education\n(Individual)",
         "Cognitive\nAbility\n(Individual)",
         "Information\nSources\n(Individual)",
         "Median Level\nof Education\n(Kecamatan)")
ndv = length(unique(df$dv))
df$modellab = rep(mods, times = rep(ndv, length(mods)))
     
df$modellab = factor(df$modellab, levels = mods)

df$category = NA
df$category[df$dv %in% c("Object: Live in Same Village", 
                          "Object: Live in Same Neighborhood", 
                          "Object: Live in Same House",
                          "Object: Worship House", 
                         "Object: Intermarriage")] = "Outgroup\nRejection"

df$category[df$dv %in% c("Prefer Muslim Candidate", 
                         "Prefer Co-Ethnic Candidate")] = "Political\nPreferences"

df$category[df$dv %in% c("Trust Fellow Muslims", 
                         "Trust Co-Ethnics", 
                         "Helping Neighbors")] = "Ingroup\nAttitudes"


df$category = factor(df$category,
                       levels = c("Outgroup\nRejection",
                                  "Political\nPreferences", 
                                  "Ingroup\nAttitudes"))


df$low95 = df$coef - qnorm(.975)*df$se
df$hi95 = df$coef + qnorm(.975)*df$se

df$col = "darkgrey"
df$col[(df$coef < 0) & (df$p<=.05)] = "red"
df$col[(df$coef > 0) & (df$p<=.05)] = "blue"

df$alph = .6
df$alph[df$p<=.05] = 1



d = subset(df, modellab == "Education\n(Individual)")
xmin = min(d$low95)
xmax = max(d$hi95)
p1 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "pt"))


d = subset(df, modellab == "Cognitive\nAbility\n(Individual)")
xmin = min(d$low95)
xmax = max(d$hi95)
p2 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == "Information\nSources\n(Individual)")
xmin = min(d$low95)
xmax = max(d$hi95)
p3 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(strip.background.y = element_blank(),
        strip.text.y = element_blank(),
        plot.margin = margin(0, 0, 0, 2, "pt")) 


d = subset(df, modellab == "Median Level\nof Education\n(Kecamatan)")
xmin = min(d$low95)
xmax = max(d$hi95)
p4 <- ggplot(d, 
             aes(x = coef, y = dv)) +
  geom_vline(xintercept = 0, lty=2) +
  geom_point(alpha = d$alph, color = d$col) +
  geom_segment(aes(x = low95, xend = hi95,
                   y = dv, yend = dv),
                   alpha = d$alph, color = d$col) + 
  scale_x_continuous(name = NULL,
                     limits = c(xmin-.02, xmax+.02),
                     breaks = round(seq(xmin, xmax, length.out = 3),2)) +
  scale_y_discrete(name = NULL,
                   labels = NULL) +
  facet_grid(rows = vars(category), cols = vars(modellab),
             scales = "free", space = "free_y") +
  theme_bw() +
  theme(plot.margin = margin(0, 0, 0, 2, "pt")) 



p <- p1 + p2 + p3 + p4 + plot_layout(nrow=1)

p
ggsave(p, filename = "Figure-8.pdf")

```


